Transport through a double quantum dot in the sequential- and co- tunneling regimes 
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We study transport through a double quantum dot, both in the sequential tunneling and cotunnel- 
ing regimes. Using a master equation approach, we find that, in the sequential tunneling regime, the 
differential conductance G as a function of the bias voltage A/i has a number of satellite peaks with 
respect to the main peak of the Coulomb blockade diamond. The position of these peaks is related 
to the interdot tunnel splitting and the singlet-triplet splitting. We find satellite peaks with both 
positive and negative values of differential conductance for realistic parameter regimes. Relating our 
theory to a microscopic (Hund-Mulliken) model for the double dot, we find a temperature regime 
for which the Hubbard ratio (=tunnel coupling over on-site Coulomb repulsion) can be extracted 
from G(Au) in the cotunneling regime. In addition, we consider a combined effect of cotunneling 
and sequential tunneling, which leads to new peaks (dips) in G(Afi) inside the Coulomb blockade 
diamond below some temperature scales, which we specify. 



I. INTRODUCTION 

In recent years, there has been great interest in few- 
electron quantum dots, or so-called artificial atoms^. 
The interest stems from a variety of fields, such as nano- 
electronics, spintronics, quantum computation, etc., 
which are all rapidly growing at present. Unlike usual 
atoms, the quantum dots can be locally manipulated by 
gates and tuned to the regimes of interest such that 
one same quantum dot can realize a whole species of 
atomic- like electronic structures (artificial atoms). A 
great deal of experimental investigation on the Kondo ef- 
fect, Coulomb blockade effect, spin blockade effect, etc., 
has been performed in recent years owing to this tunabil- 
ity feature. Yet not a less important feature of quantum 
dots is the possibility to architecture and control, to a 
great extent, the coupling to their surrounding, be it a 
dissipative environment, classical (gate) fields, or other 
quantum dots. This makes quantum dots attractive for 
quantum computation. 

A most promising candidate for qubit (quantum bit) 
in solid state is the electron spin 5 . Confining electrons to 
quantum dots naturally defines the qubit as the spin up 
and spin down states of the dot, provided the dot contains 
an odd number of electrons. At the ultimate level of 
control over the electron charge, the quantum dot can be 
tuned to confine one single electron^, thus implementing 
the artificial version of the hydrogen atom. This has 
been successfully achieved in recent years, first in vertical 
dots^ and lately also in lateral dots4, due to a special 
design of top gates. Furthermore, observation of shell 
filling of dot orbitals and the Hund's rule in symmetric 
dots^ indicates that the electron spin is a well defined 
and relevant degree of freedom in few-electron quantum 
dots and that achieving control over it should be feasible 
experimentally in the near future. 

A very recent step forward in accessing the electron 
spin in quantum dots was made in Refs.0 andllCt where 
the Zeeman splitting of an electron in a lateral quantum 
dot has been measured directly by means of dc transport 



spectroscopy. Resolution of the Zeeman sublevels in a 
magnetic field B > 5 T allowed the authors of Ref . to 
measure the spin relaxation time T\ for a single electron 
by means of a pulsed relaxation measurement technique. 
This techniqueii uses a sequence of pulses, which allow 
for sequential-tunneling into excited states of the quan- 
tum dot 12 . Relaxation from such an excited state can be 
monitored in the average current versus the pulse tim- 
ing, see Refs. Eland© for details. Observation of an or- 
bital relaxation time of 10 nsii and a Ti-time larger than 
50 /is at B = 7.5 T 9 , illustrates the long-lived nature of 
the electron spini£. 

Integrating qubits into a quantum computer is at 
present at the stage of realizing a two-qubit circuit. In 
the language of artificial atoms this translates to ma- 
nipulating a hydrogen molecule. Here, the spins on two 
dots interact with the Heisenberg exchange interaction 
i?gp m = JSi • S2, where Si^ are spin 1/2 operators, and 
J is the exchange coupling constant. Achieving control 
over the coupling constant J is as important for the qubit 
as controlling its Zeeman energy for spin quantization 
along two independent axes. The latter allows for arbi- 
trary single-spin rotations on the Bloch sphere, whereas 
a combination of the exchange interaction and single- 
spin rotations 5 " allows one to create non-local quantum- 
mechanical correlations between the qubits, often re- 
ferred to as entanglement. Moreover, as demonstrated in 
Ref.H the exchange interaction together with single-spin 
rotations suffices for universal quantum computation. 

The challenge of implementing an artificial Ho 
molecule has been met in recent works, see Refs. hj 
and El In typical structures, transport measurements 
alone do not always give reliable information about the 
number of electrons on each dot. The current through 
the double dot (DD) is often too small to be measured 
and the tunnel contacts are likely to be pinched off in the 
few-electron regime, due to a high level of depletion in 
the dot region caused by the electrostatic potential of the 
gates. A complementary method of charge control was 
used in Ref. which came with integrating the quan- 
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turn dots with sensitive charge detectorsiSiii — quantum 
point contacts (QPCs) placed nearby the quantum dots 
and tuned to half pinch off. Such a charge detector is 
capable of sensing a change in the dot occupation num- 
ber by as much as a fraction of an electronic. Using this 
method of charge control should allow one to concentrate 
in more detail on the transport properties of the DD. As 
we show in this paper, the differential conductance of 
the DD as a function of the bias voltage provides valu- 
able information about the DD parameters, including the 
exchange constant J. 

We begin with a Hund-Mulliken model for the DD in 
Section [H] We introduce the parameters which char- 
acterize the DD. Next, throughout the paper, we refer 
to these parameters as phenomenological ones, seeking 
ways to extract them from possible experiments. In Sec- 
tion ^0] we use a master equation approach and derive 
the sequential-tunneling current 7 at a finite bias voltage 
Afi. Here, we find a novel feature for the DD: the differ- 
ential conductance G = edl /dAfi as a function of A/i has 
a peak of negative G for typical DD parameters. We find 
a number of additional peaks, which all together allow 
one to extract the exchange constant J and the DD tun- 
nel splitting 2to. In Section llVl we consider cotunneling 
through a Coulomb blockaded DD with one electron. We 
calculate the elastic and inelastic components of the cur- 
rent and specify a temperature regime of strong heating, 
where the asymmetry of the coupling to the leads, 77, can 
be extracted from transport measurements. The proce- 
dure is explained in detail in Section^ where we consider 
the cotunneling through the Coulomb blockaded DD with 
two electrons. We show that the interaction parameter (ft, 
which gives the Hubbard ratio and entanglement between 
the electrons in the DD singlet state, can be extracted 
from transport measurements. Finally, in Section IVTl we 
consider a combined effect of sequential tunneling and 
cotunneling for the DD. Here, again we find additional 
peaks or dips, which can occur in the Coulomb blockaded 
valleys for the non-linear conductance below certain tem- 
perature scales, which we specify. 



II. ENERGY SPECTRUM OF A DOUBLE 
QUANTUM DOT 

For definitcness we consider lateral quantum dots, 
which are usually formed by gating a two-dimensional 
electron gas (2DEG) under the surface of a substrate. 
The 2DEG is depleted in the regions under the gates and, 
with an appropriate gate design 4 , one can achieve a de- 
population of the dots down to 1 and electrons per dot, 
avoiding pinching of the dots from the rest of the 2DEG. 
The low energy sector of a DD at occupation number 
N = Ni + N2 = 1 consists of two tunnel-split energy 
levels, which we label by the orbital quantum number 
n, with n — + standing for the symmetric orbital and 
n = — for the anti-symmetric one. The states of the DD 



with one electron can then be written as 



= 4 jo) 



(i) 



where a denotes the spin degeneracy of each level (we 
neglect the Zeeman splitting) , d^ a is the electron creation 
operator, and |0) is the DD state with zero electrons. The 
splitting between the two levels is given by 2to, where to 
is the interdot tunneling amplitude. We assume weak 
tunnel-coupling between the dots such that to <C ftiuo, 
where ftwo is the size-quantization energy of a single dot. 
Then, for the occupation number N = Ni + N 2 = 2, the 
lowest energy states, one singlet state and 3 triplet states 
are given by 
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\T Q ) = -L(d^dl l +dl l d^)\0). 



(2) 



The splitting between the singlet and the triplet, J — 
E\t) — E\s)i pl a ys the role of the Heisenberg exchange 
interaction for the two electron spins in the DD, H s ^ m = 
JSi • S2. The interaction parameter cj>, entering the sin- 
glet state in J2J , is determined by a competition between 
tunneling and Coulomb interaction in the DD, and it can 
be calculated 18 (Hund-Mulliken method) to be 



(3) 



where tn and Uh are the so called extended inter-dot 
tunneling amplitude and on-site Coulomb repulsion, re- 
spectively*^. We note that tn = to + tc — to, where tc is 
a Coulomb contribution, which vanishes with detaching 
the dots (to — * 0). To illustrate the meaning of cf>, we 
also present the states in terms of orbitals, which are 
mostly localized on one of the dots, see Appendix A. The 
double occupancy in the singlet state is given by 
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2(1 + ^2) 



(4) 



The parameter (j> also determines the entanglement be- 
tween the two electrons in the singlet state. While <f> can 
be used on its own as a measure of entanglement, we are 
presenting here a formula for the concurrence^ of the 
singlet \S), 
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(5) 



The entanglement in the \Tq) state is maximal (c = 1) at 
all values of to. 

The energies of the considered DD states are given by 

(6) 
(7) 

E {s) =2i+U 12 - J, (8) 
E\t) = 2e + U 12 , (9) 



^|o> — 0, 

E\ n <j) = £ - nt-a, 
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where e ~ hu>o — V g and U12 — e 2 /2Ci2, with V g being 
the energy shift due to the common gate potential V g 
(see Fig.nj and C%2 the mutual capacitance between the 
dots. 

For small dots (ftujo > Uh), the next excited states are 
the following two singlets 



1 



(10) 



|52) 



with Ei si) ~ Ei 



- 7 ===(^4 T 4 i + dL T di i )io), (11) 



\S2) '- 



E\s) + Uh- The states (01 together 



with 11 0|l and (|llfl complete the resolution of unity for 
two electrons in the DD orbitals n = ±. 

Finally, we note that according to the Hund-Mullikcn 
method the exchange constant J consists of two compo- 
nents, 



J = Vc + Jh, 



(12) 



where Vc < is responsible for a singlet-triplet transition 
at a finite magnetic field, see Ref. fl9l and 



Jh 



1 



1, 



, v U% + 16^ - \U H 



(13) 



resembles the exchange constant obtained in the stan- 
dard Hubbard model for on-site Coulomb repulsion. For 
weakly coupled quantum dots, we have Jh ~ ^h /Uh- 



III. SEQUENTIAL TUNNELING THROUGH 
THE DOUBLE DOT 

The setup we are considering is shown in Fig. ^t. Each 
dot is tunnel-coupled to a metallic lead via a point con- 
tact, forming a series lead-dot-dot-lead setup [for a paral- 
lel configuration see R.efs. I21I and l22l| . The point contact 
in lateral structures has a smooth (wave-guide-like) po- 
tential, providing a number of lead modes (channels) that 
can couple to the dot. By constraining the point contact, 
the lead modes can be pinched off one by one, due to the 
transverse quantization in the contact region. When the 
last mode in each point contact is about to be pinched 
off, the structure shows Coulomb blockade (CB) effect at 
low temperatures. In this regime, the point contacts can 
be treated as tunnel junctions, with tunneling amplitudes 
and tfi, for the left and right dot, respectively. The 
leads are then single-channel Fermi-liquid leads, which 
we describe by the Hamiltonian 



leads 



(14) 



I ha 



where c lka creates an electron with momentum k (en- 
ergy £fc) and spin a in lead I = L, R. The tunneling 
between the DD and the leads is described by the tun- 
neling Hamiltonian 



H, 



J2 H t 

1 



l nka 



h.C. 



(15) 
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FIG. 1: a) Two coupled quantum dots, with the inter-dot 
tunneling amplitude to, attached to metallic leads at differ- 
ent chemical potentials (il and [ir = /il — A^i. The tun- 
nel coupling between the dots and leads is characterized at 
to — by the tunneling amplitudes tL and ta. b) A schematic 
plot of the linear (A^i —> 0) conductance G through the dou- 
ble dot (DD) as a function of the gate voltage V g , showing 
the Coulomb blockade (CB) effect, with sequential-tunneling 
peaks and CB valleys. The number of electrons contained 
in the DD, N = Ni + N2, is fixed in the valleys between the 
peaks. The position away from the TV = 1, 2 peak is measured 
by 5E, which takes on negative (positive) values in the N = 1 
(TV = 2) CB valley. 



with the tunneling amplitudes 
t L 



a/2(1±5) 



t R ,± = ± 



v/2(l±5) 



(16) 



Here, S ~ to/hujQ is the overlap integral between the two 
dots orbital wave functions (S — (<Pl\<Pr))- Formulas 
l|lfi|l account for the hybridization of the DD orbitals at 
2io ^ fi^o with the accuracy of the method of molecular 
orbitals (Hund-Mulliken method). The tunnel- coupling 
to the leads broadens the DD levels n = ±, introducing 
the level width T n = irv(\tL n \ 2 + |<_r„| 2 ), where v is the 
lead density of states per spin. For our convenience, we 
also use the notations: Ti — irv\ti\ 2 : with I — L,R, and 
V = K-rIVI^lI 2 - Throughout the paper, we assume r± < 
T, where T is the temperature. 

The usual CB stability diagramiS for a DD repre- 
sents a honeycomb structure of increased linear conduc- 
tance plotted versus V g i, V g i — the gate voltages con- 
trolling each of the two dots, respectively. We are in- 
terested in the case when the two dots are similar, and 
therefore consider the diagonal of the stability diagram, 
V g i = V g i = V g , in the vicinity of (1,1)-(1,0)-(0,1) triple 
point. However, we will be interested in large applied 
bias voltages A/i = — /j,r, assuming that the bias 
voltage drop occurs on the structure as a whole and does 
not shift it away from the diagonal of the stability dia- 
gram. The bias voltage A/i can be applied in different 
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ways, e.g. equally distributed between the left and right 
leads, or applied to one of the leads alone. To cover all 
possibilities, we assume the chemical potentials of the left 
and right leads to be respectively equal to: 



It is convenient for the calculation to trace out the 
spin degeneracy of and of the triplets ©■ We map 
the degenerate levels onto non-degenerate ones, using the 
following replacement^ 



p L = p + Ap L , 
fi R = p, - Ap R , 



(17) 



1 

Ni 



iei 



fi 



FI- 



(22) 



with ApL+Ap R = Ap, and the bias-asymmetry parame- 
ter: a = Ap R / ApL,. The position in the CB is controlled 
by the common gate voltage V g . In the linear regime 
(Ap — > 0) the conductance G as a function of V g shows 
peaks at the degeneracy points of the chemical potential 
in the DD and in the leads, see Fig. We focus on the 
sequential tunneling peak, where the number of electrons 
in the DD fluctuates between N = 1 and N = 2. The de- 
generacy condition is given by: E(2) — E(l) = p, where 
E(N) is the ground state energy of the DD with N elec- 
trons. We assume the singlet (jHJ to be the ground state 
for TV = 2; therefore, £7(2) = E ]S ) and £7(1) = E\ +i<T ). 
The position in a CB valley is characterized by the addi- 
tion/extraction energy: E± (N) = E(N ± 1) - E(N) + p. 
The distance away from the peak (in the scale of V g oc V g ) 
is measured by 



where Nj is the degeneracy of level /. Thus, from here 
on, we deal with 4 non-degenerate levels, denoted as |+), 
|— ), \S), \T). The transition rates between these states, 
for I = L, are: 



2p^l f{ SE + A, L) , 



wt 



•f(—SE — 2to — A/xl), 



1 

2^^^f(5E + 2t Q + Ap L ) 
3 2tt 



SE = E\ 



+,<?) 



E, 



s) + M- 



(18) 



h -v\t L ,-\ J f(J-6E-Aii L ), 

W +' T = Y^ tL ^ 2f ^ J + 5E + Aml) ' 



Positive (negative) values of SE correspond to positions 
in the N = 2 (N = 1) CB valley. The size of the N = 1 
CB valley (distance between its 1st and 2nd CB peaks) 
is given by E + (l) + £_(1) = U 12 + 2t - J. The size of 
the N = 2 CB valley is of order of U H (U H >U 12 ). 

The sequential tunneling through the DD is described 
by the golden rule rates: 



Wt. 



~v\t L , + \ 2 f(J-SE-2t -Ap L ), 



2tt 

Wt, T = —„\t Li+ \ 2 f(-J + 5E + 2t + Ap L ).(23) 



W l Mm = ^-T l \( M ' n \ C lk,H. 



Expressions for / = R are obtained from above by sub- 
stituting L — > R and Ap^ — > — A fin. Here, f(E) = 
1/ [1 + exp(E/T)] is the Fermi function. 



n.ka 



8{Sk ~ EMm)pf. 



W l mM = ^ K m; n \cikaH l T \M ; n) \* 8(e k + E Mm )pf. 



(19) 
(20) 



A. Master Equation 



where m stands for one of the states Q , and M for one 
of the states © ; W l Mm is the probability rate for the DD 
transition from m to M by exchanging an electron with 
the lead I; Eum = Em — E m with Em, E m being one of 
© — ©■ The averaging in fl§|) and (|2T))) is performed 
over the Fermi-sea states |n) with the grand canon- 
ical density matrix pf = Z^ 1 cxp(— A^/fceT), where 

Zi = Tr ; exp(-Ki/k B T), K t = H- pi ^ fc(T cj kcr ci k a, and 
T is the temperature (we set k R = 1 in what follows). In 
(|19|l and l]2Up. we use the notations: pf n = (n\pf\n) and 
|m;n) = |m)|n). As a consequence of the thermal equi- 
librium in the leads, the rates l|19|) and H2U|) are related 
by the detailed-balance formula: 



Assuming a large temperature T 3> r± , but still much 
smaller than the scales of interest, we describe the DD 
state by the probability p p for the DD to be in the level 
pG {|+), |— ), \S), \T)}. In the stationary limit, p p obeys 
the balance equations: 

(W s ,+ + W T ,+) P+ = W+^sPs + W+^tPt , 
(Ws- + Wt,-) p- = W-^sPs + W^^tPt , 
(W +jS + W-, s ) ps = W s ,+p+ + W s ,-p- , 
(W +;T + W_ t) Pt = W T ,+p+ + W T ,-P- , (24) 

where we used the notation W pp i = W pp , + W pp , . Only 
three of the equations H24J) are linearly independent. 
Choosing any three of them, and using the normaliza- 
tion condition 



W} nM = W l Mm exp[(E 



Mr 



(21) 



P+ 



Ps + Pt = 1 



(25) 



We provide explicit expressions for the rates W \ 
Appendix B. 
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one can find the solution for p p . However, it is convenient 
for the further discussion to describe the non-equilibrium 



in the DD by the following balance ratios: 

Ps + Pt 

T = 1 ' 

P+ + P- 

P = Pt/ps , 
7 = P-/P+ , 



(26) 

(27) 
(28) 



which give the occupation probability p p of the DD states 
as 



1 



P+ = 



Ps = 



(l + r)(l + 7 ) ' 



(l + r)(l+/3) ' 



P- = 



Pt = 



7 



(l + r)(l + 7 ) 
r_l 

(l + r)(l+j9) 



(29) 



Expressions for r, (3, 7 are given in Appendix C. In the 
linear regime, the DD is in thermodynamic equilibrium, 
and the occupation of the states is determined by the 
temperature T. For this regime, we find the equilibrium 
values: 
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FIG. 2: (a) The ratio r = (ps + pt) / (p+ + P-) plotted versus 
Ap/2 = ApL,R for SE < 0. (b) The inverse ratio 1/r versus 
Ap/2 for SE > 0. For both figures, we choose the units of 
energy such that \SE\ — 1. We use the following parameters: 
2t = 0.8, J = 0.2, S = 0.5, and = 0.4. The grids above the 
figures (a) and (6) show the division into the intervals 1341 
and 141H . resp. The insets to the figures show the transitions 
which switch on at Ap/2 = \SE\. 



(3 T = 3 exp(- J/T), 7 T = exp(-2t /T), (30) 
t (l + f3 T )exp(SE/T) 

T = WTln ' (31) 

In the non-linear regime, the deviation from these equi- 
librium values due to the applied bias describe the heating 
effect in the DD. To simplify our further considerations 
we make the following assumptions: 

a) \SE\ >2t > J>0, (32) 

b) T< J, 2t - J, \SE\-2tQ. (33) 

For dcfinitcness, we also assume a symmetric bias situa- 
tion with A/i£ = ApR = Ap/2. Next, we consider the 
two following cases. 

1. Sequential tunneling on the N — 1 CB valley side 
(SE<0). 

As mentioned in Section [HI the N = 1 CB valley has 
the size U12 + 2t n — J . We can estimate t/12 for a lateral 
structure from the Hund-Mulliken method in the absence 
of screening of Coulomb interaction due to top gates. We 
obtain^ U12 ~ 2.7 meV for GaAs quantum dots with 
Hujq = 3 meV coupled so that the distance between the 
dots centers is ~ 2a b, with as = y/h/muiQ being the 
Bohr radius of one dot. For the on-site Coulomb re- 
pulsion we obtain Uh ~ 4.5 meV. The screening from 
the top gates, which depends on the design of the struc- 
ture and on the thickness of the insulating layer be- 
tween the 2DEG and top gates, reduces, in practice, 
the inter-dot Coulomb repulsion as compared to the on- 
site one. However, we still assume a sizable U12 such 
that U12 > 2to, and thus we can neglect the contri- 
bution from the N = 0, 1 sequential-tunneling peak at 
An/2 ~ \SE\ + J < (t/12 + 2t - J)/2. 



For further consideration it is convenient to divide the 
range of the applied bias into the following intervals: 

I. < Ap/2 < \SE\ - 2t , 

II. \SE\ ~ 2t < Ap/2 < \SE\ - 2t a + J , 

III. \5E\ -2t + J< Ap/2 < \8E\ , 

IV. \SE\ < Ap/2 < \5E\ + J , 

V. \SE\ + J < Ap/2. (34) 

At temperatures satisfying l|33l) , the solution of the mas- 
ter equation l|24|l is constant within each of these inter- 
vals. We plot the quantities (|2*rj|l . (|2*7|l and 1(25)1 versus 
Ap/2 in Figs. and^Ji, respectively. The units of 

energy are chosen such that \8E\ — 1. Fig. [2Ji shows 
the balance between the DD being in the sector with 
N = 2 electrons and the DD being in the sector with 
N = 1 electrons. At small bias voltages (I, II and III), 
the occupation of the N = 2 sector is suppressed as 
exp [— (\SE\ — Ap/2)/T], because of the Coulomb block- 
ade in the N = 1 valley. At the bias voltage Ap/2 = \SE\, 
the left lead chemical potential p^ reaches the value of 
the DD N = 2 chemical potential. At this point the 
following sequence of transitions becomes possible: 

|+) <=* |S) «=* I-) ^ |T), (35) 

which changes the population probabilities p p in the DD. 
It is important to note the difference between the inter- 
vals IV and V. In the interval V we get universal results: 
t = 1/rj, f3 = 3, and 7 = 1. This corresponds to set- 
ting T — > 00 in Eq. (|30l) for [3 and 7, and recognizing 
that (ps + Pt) ~ and (p + + pJ) ~ Tr, which yields 
t = Tl/Tr = 1/r] according to Eq. In contrast, in 

the interval IV the plateau values of r, j3 and 7 depend 
on the DD parameters. For example, in Fig.[2Ji this non- 
universal value of t is denoted by Tjy, and we find that 
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FIG. 3: With the same parameters as in Fig. the ratio 
P = Pt/ps versus Ap/2 for: (a) 5-E < and (6) SE > 0. 
The shaded area in (a) shows the inactive region with no 
population of the N = 2 sector, see Fig. [5Jj. We note that 
this region can be active if a combined effect of sequential 
tunneling and cotunneling is considered, see Sec. lVIl 



yr< = r, 



1 + l/> 2 + 2 X (2 + 2 ) 
l + l/0 2 + 2x73 + 877/3(1 + 7?) 



. (36) 



where x = T-/r+ = (1 + 5)/(l - S). This "universal- 
ity versus non-universality" depends on the way the se- 
quence (|35|l is closed. For the interval IV, only the transi- 
tion |T) — * |+) is allowed, whereas the reverse transition 
is forbidden by energy conservation. For the interval V, 
however, the sequence is closed by \T) ^ |+). We can 
express the fact that, for the regime on the plateau V, the 
results are universal by formulating a principle of detailed 
balance. For the rates entering the master equation (|24H . 
such a non- equilibrium detailed-balance principle can be 
written as follows 



W s .± = tt w ±,s, 
2rj 



2,/ 



W± : 



(37) 
(38) 



This suffices to obtain the universal result for r, (3 and 7 
from the master equation H24(l . 

For the Fig. EJi, we find that /3 changes from zero 
to a value /3 IV at the border between the intervals II 
and III. This is related to the fact that the transition 
\S) — > | — ) — > \T) could occur, provided there was a non- 
vanishing population in the N — 2 sector. But since 
t = for this interval, see Fig. |2Ji, this step-like change 
in (3 will not be observed in physical quantities. Interest- 
ingly, when the N = 2 sector acquires non-zero popula- 
tion at the border between the intervals III and IV, the 
ratio (3 = Pt/ps stays constant as a function of Ap, and 
T, for T satisfying l|33|l ■ The value of f3^ is given by 



i +ci + ,: x -|(i+.i/^ 2 ) 



(39) 



FIG. 4: With the same parameters as in Fig. [5] the ratio 
7 = p-/p+ versus Ap/2 for: (a) SE < and (b) 5E > 0. The 
inset in (b) shows the sequence of transition which takes place 
at Ap/2 > 5E + 27jo — J (compare with the inset in Fig.|5[)). 



and that of 7^, by 

1 / < -1,3 l + l/0 2 + x(l + 2 ) 

1/7lv ^ 1+ 2 x + n/ii + r,) ■ (40) 

The population probability of each of the levels can be 
obtained from the formulas H29|l. We plot p + , p_, pg, 
and pt on Fig. [SJi for the above discussed situation. 



Sequential tunneling on the N 
{SE>0). 



2 CB valley side 



Here, the relevant intervals of applied bias are: 

I. < A/j/2 < SE- J, 

II. SE- J < Ap/2 < SE, 

III. SE < Ap/2 <SE + 2t a - J , 

IV. SE + 2t — J < Ap/2 < SE + 2t Q , 

V. SE + 2t < Ap/2 , (41) 

and we assume 8E + 2to < Uh/2 such that the DD is not 
populated with 3 electrons while raising pl ■ We plot 1/r , 
f3 and 7 versus Ap/2 in Figs.[2),[It and^J), respectively. 
At Ap = the DD is in the iV = 2 CB valley and as in 
the previous case a sizable change in the DD population 
occurs when Ap/2 = SE. At this point the chemical 
potential pr in the right lead is low enough such that an 
electron from the DD with N = 2 can occupy an empty 
place above the Fermi sea in the right lead. The following 
sequence of transition is immediately activated 



I S) 



(42) 



which yields universal values: r = 2/77, /3 = 3 and 7 = 0, 
in the interval III. The corresponding principle of detailed 
balance is obtained if we disregard the level |— ) in the 
master equation (|24|l and use l|37(l and (|38J) with ± — > 
+. The left-to-right processes of the sequence lf4*2*|) arc 
illustrated in the inset to Fig. Yet two other changes 




FIG. 5: The population probabilities p p versus Ap/2 calcu- 
lated from Eqs. 121 H with the values of r, j3 and 7 given by 
respectively Figs^EJandglfor: (a) SE = -1 and (6) SE = 1. 



FIG. 6: The bias dependence of the differential conductance 
G for: (a) SE = — 1 and (b) 5E — 1. The ordinate axis is 
scaled by Gi = e 2 irv\ti\ 2 /hT . We use the same parameters as 
in Fig. |5|and with t L = t R (G L = Gr). 



in 1/t occur at Ap = SE + 2t — J and Ap = SE + 2t , 
see Fig. [2>. The value of 1/t on the plateau IV in Fig.[3> 
is given by 



'IV 



X (7 + 30 2 )/4-l/2 



3(l + l/0 2 )/(l + 77) + l + X 



(43) 



Fig. 13* shows that, in the interval IV, the triplet level 
loses its population probability relative to the singlet 
level, with the value of (3 being given by 



i/Av = ~ 



1 



2 + 3 X (l + 2 ) 



2 X + 3(l + l/0 2 )/(l + r ? ) 



(44) 



Just similarly to Section lill A II we have here the follow- 
ing sequence of transitions 



\S) 



\T) 



(45) 



which is closed in the interval IV by |— ) — > IS"), and 
in the interval V by |— ) T± \S). The latter results in 
t = 1/ry, /3 = 3, and 7 = 1, which is identical with the 
universal result in Section fill A II The detailed balance 
for this case is also given by Eqs. (|37|l and (|38|l . Finally, 
the non- universal value of 7 in Fig. ^ is given by 



V7fv = l+o 



2 (l + r,)(l+ X ) 



1 + 1/4> 2 



(46) 



and we present each of the population probabilities p p 
for the considered case in Fig. [3J>. 



B. Sequential Tunneling Current 

The electron (particle) current, flowing from the DD 
into the lead I, reads 



I 1 = (W[, s + wL tS ) p s + (w| iT + Wi tT ) PT 



(w l St+ + w l Tt , 



w 



(47) 



In the stationary regime described by l|24|l , one has I L — 
—I R = I/\e\. The differential conductance G = edl/dAp 



as a function of SE and Ap can be evaluated for different 
regimes of interest. 

For the regime studied in Section UlI A II we can use a 
simplified formula for the current, namely 



I/Io R 



1 + (j) 2 v 1 + s i-s 



PS 



1-S 2 



p T , (48) 



where = \e\T R /H. We plot G vs Afi/2 for SE = 
-1 in Fig. Eh The main peak at Afi/2 = \8E\ = 1 
acquires a satellite peak at A/i/2 = \SE\ + J, which for 
the given parameters has a larger amplitude than the 
main peak. The origin of the satellite peak is closely 
related to the heating effects discussed in Section llll A II 
Eq. I|48|) shows that the changes in ps and pr as functions 
of A/i are directly reflected in the current /. 

For the regime studied in Section llll A 21 we can use a 
simplified formula for the current, namely 



I/Io L 



SI 
1 



1-S1 



3/2 
1-5 
3/2 



p- 



(49) 



where J L = \e\T L /h. We plot G vs A/i/2 for SE = 1 
in Fig. [U>. The main peak at Ap,/2 = SE = 1 ac- 
quires two satellite peaks at Ap/2 = \SE\ + 2to — J and 
A/z/2 = |<5_E| + 2to. Interestingly, the first satellite peak 
has negative differential conductance for the given param- 
eter values. Eq. I|49l) shows that the current / reflects the 
changes in p + and p_ as functions of Ap, discussed in 
Section IIII A 21 The negative value of G is due to the 
decrease of p+ when going from the interval III to the 
interval IV (see Fig. 03>) and different tunnel coupling to 
the n — + and n = — energy levels. At the very origin of 
negative G lies the Coulomb interaction in the DD, which 
allows us to consider a truncated Hilbert space, namely, 
consisting of the states Q and (0) . 

Using Eqs. (|47|l and l|29|l . we calculate the differential 
conductance G for the whole range of variables SE and 
Ap. Fig. [7| shows a gray-scale plot of G for the case 
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0.5 






5E 



FIG. 7: A gray scale plot of G versus SE and A/i/2 = ApL = 
ApR. The white (black) color corresponds to positive (neg- 
ative) values of G; gray stands for G = 0. Here, we use: 
2t = 1, J = 0.25, 5 = 0.6, = 0.3 and r/ = 1. 



of symmetric biasing: A/x^ = A/i^ = A/z/2. The gray 
color corresponds to G = 0, the white (black) color cor- 
responds to positive (negative) values of G. We note 
that the black line on Fig. [7| terminates at a satellite line 
(5E = —to + J, Ap/2 = to), unlike the other two white 
(satellite) lines, which terminate at the main sequential- 
tunneling peaks. This can be attributed to the origin 
of the black line: change in the rate for excited state to 
excited state transition, see \T) — > |— ) in the inset of 
Fig.© (cf. inset of Fig.©). For A^/2 < t , the tran- 
sition |T) — > |— ) is blocked due to energy conservation. 
We note that "excited state to excited state" sequential- 
tunneling satellite lines have been observed experimen- 
tally for single dotsi£. Finally, in Fig. |SJ we present a 
gray-scale plot of G for the case of asymmetric biasing: 
AfiL = Afi and Ap^ = 0. 



C. Charge Detection via QPCs 

Using quantum point contacts (QPCs) placed in the 
neighborhood of the quantum dots provides an additional 
channel of information about the DD 1 ^. Here, we con- 
sider the average charge on the DD, (N) = p + + p- + 
2(ps + Pt)- With the help of Eqs. and w e relate 
(N) to the parameter r as follows 



(TV) 



1 



1 



(50) 



At a large bias voltage, corresponding to, e.g., the in- 
terval V in Section IIII A II see Eq. J33J, the DD oc- 
cupation number fluctuates between 1 and 2, being on 
average (N) = 1 + 1/(1 + i]), for Ap > 0. This rela- 
tion can be used to find out the asymmetry parameter 
i]. In the case of symmetric coupling to the leads one 
has (N) = 1.5. In the interval IV of Section IhTXTI 
the average (N) assumes a non-universal value, deter- 






8E 



FIG. 8: The same as in Fig|7| but for asymmetric biasing 
Ap = ApL, Ap R = 0. 



mined by Eq. I|36|l . This result can, in principle, be 
used to extract the overlap integral S [or equivalcntly 
X = r_/r+ = (1 + 5)/(l - S)] in the case when the 
parameter <f> is known, e.g. as from Section^ For sym- 
metric coupling to the leads (rj = 1) and in the limit of 
weakly coupled dots (0 — ► 1, S — > 0), Eq. l!3rj|l yields 
t = 1/3, or equivalcntly (N) ~ 1.25. This can be dis- 
tinguished from the value (N) = 1.5 in the interval V, 
with the available sensitivity of QPCs, which is ~ 0.1 e, 
see Ref. Q. Thus, we conclude that the satellite peaks 
obtained for the differential conductance in Section lill Bl 
(Figs.0and[Hl can be observed also using QPCs, 24 Next, 
we discuss the back action of the QPC on the DD. 

The QPC detects the charge of the quantum dot, that 
it is attached to, via the electrostatic potential the dot 
induces in the QPC region. We consider the DD and one 
QPC close to the right dot. The term responsible for 
relaxation of the DD states is as follows^ 



SHqpc = -ife(i)V4.<T d -»,<r 



(51) 



where Se(t) is the fluctuating field of the QPC, with 
Se(t) — (bar denotes average over the QPC degrees 
of freedom) . The relaxation rate is given by2& 



1 



J(E if /h), (52) 



where £Jj/ = Ei — Ef is the transition energy. It follows 
from Eq. I|52|l that the QPC can induce relaxation only 
between the states |+, a) and | — , a); all other transitions 
are forbidden due to the spin and charge conservation 
in the Hamiltonian l|51() . The spectral function J(uj) is 
defined as follows 



J(uj) = / Se(0)Se{t)e- lut dt. 



(53) 



The QPC field, Se(t) — e(t) — e(t), is proportional to the 
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charge density at the QPC, 



s(t) = 



IVkk'v 



^LWi'k'At), (54) 



where the indices I, I' — C, 1Z denote the two leads of the 
QPC (not to be confused with the DD leads), fn is the 
chemical potential of lead I, flc = fin + A/j, and cj kcr (t) 
creates a QPC electron. For u>, Afi/h < ui c , we have 



the DD in a state with one extra (fewer) electron and is 
transferred to either lead, leaving the DD with the same 
energy (elastic cotunneling) or with a different energy 
(inelastic cotunneling). The cotunneling rate for a lead 
electron to go from lead I to lead V and the DD to go 
from state m to state n is given by the golden rule rate^ 

w nm (l',l) = ^^|(n|T r; |m)| 2 <5(£ mn + A M „0P^ (59) 



J(w) 



= Airhv 2 (X 2 CC + \ 2 nn ) 
+Anhv 2 X 2 cn [<d(huj 4 



e(fiw) 

A/5) + Q{hw 



Am)], (55) 



where Q(E) = Ej (1 - exp(-E/T)), v is the density of 
states per spin in the QPC, and uj c is the high-frequency 
cutoff (order of bandwidth). Formula 1(52(1 is valid for 
weak coupling, i.e. uXw -C 1. At fku,T <C A/2, Eq. i(5B|) 
reduces to J(0) = ^irTuPX 2 ^ |A/2|, which formally coin- 
cides with J(0) = HTd in the weak coupling limit, where 
is the QPC decoherence rate for a single do1*2L2£. We 
summarize here by writing down the non-zero relaxation 
rates due to the QPC: 



W+- = ±J{2to/h), 
W„+ = ±.J(-2to/h), 



(56) 
(57) 



where we have already summed over the spin degeneracy, 
using the rule We note that, for \Ap.\ < 2t , the 

rate VF_ + is exponentially suppressed at low tempera- 
tures. The rates (|56|l and l|57|) describe relaxation of the 
DD due to, respectively, excitation and annihilation of 
an electron-hole pair in the QPC. 

Including the rate W-\ into the balance equations 124(1 

yields a correction, 5N, to the average occupation num- 
ber on the DD, (N') = (N) + 5N, where (N) is given by 
Eq. 1(50(1 . In leading order in W-\ , the correction reads 



8N 



1 



7 



W+ 



(1 + r) 2 (l + 7 ) 2 W 



(58) 



where 1/W is given by Eq. I(C4(I. Note that SN is pro- 
portional to 7 and, therefore, vanishes in the interval III 
of Section IlII A 21 see Fig.|3J>. In contrast, SN is finite in 
the interval V of the same section. Thus, if the average 
N is measured in the interval III, its value is given by 
(N') = 1 + 1/(1 + 77/2), see Fig.© and Eq.[5U| whereas, 
if the average N is measured in the interval V, its value 
is given by (N') = 1 + 1/(1 + rf) + SN. 



IV. COTUNNELING IN THE N = 1 COULOMB 
BLOCKADE VALLEY 



The CB valley with N = 1 has the width £L(1) + 
E + (l) = U12 + 2t — J- Cotunneling dominates the con- 
ductance in the valley at T <C U12 + 2t — J. In a cotun- 
neling process, a lead electron (hole) coherently occupies 



where A/in' = fJ-i — Pv and E mn = E m — E n . Here, we 
use the notation E m = E m + E 7fl and |m) = |m)|m), 
where E m and \m) denote the DD energy JJJ and state 
||TJ, resp.; furthermore, |m) is an eigenstate of Kl + Kr 
with energy E r n . The averaging in 1(59(1 is performed over 
the leads at thermal equilibrium with the density matrix 
= ® pji\fh). The stationary (u — > 0) cotun- 

neling is described by the effective T-matrix amplitudes 
(second order of perturbation theory) 



E.* . I d} lcr d n i a > j 



+ E 4k>a> C lk* J ,(60) 

where we assumed J <C E+(l) <C Ujj and used A = 
|jS}(i?| + J2i I^K^I to exclude virtual transitions to the 
states (HHJl, {HI- Formula is valid for A/t,T < 
E±(l). Similarly, to Sec. IIIII we trace out the spin de- 
generacy in (|TJ, using the rule 1(22(1 . The cotunneling 
rates can then be presented as follows 



>(J'0 = ^ 2 e(E n 

n 



A m >)M^ m , 



(61) 



where Q(E) = E/ (1 — exp(— E/T)) and M l n l rn are given 
in Appendix [D] The state of the DD is described by 
p+ = 1/(1 + 7) and p_ =7/(1 + 7), with 



7 = 



6(-2to + A/jl) + O{~2t - A/i) + fa + l/7 ? )e(-2t ) 
9(2i + Afi) + 9(2i - A M ) + (ry + l/7?)6(2* ) : 



(62) 



where ui„ m = w nm (l', I). We note that this result is 
universal and does not depend on the number of virtual 
states taken into account. The cotunneling current is 
given by 



e w nmPr> 



(63) 



where w„ m = w nm (R, L) - w nrn (L, R). For T < 2t , we 
define the elastic and inelastic components of the current, 
I = I c i + iinci, as follows 



"++> 



(64) 



/inci = ewL + p+ + e (w I + _ + w I __ - w J ++ ) p_, (65) 



where we used that ~}2 n Pn = 1- We note that the com- 
ponent / c i is a linear function of A/i, whereas dL mc \/dAfi 
has a step-like A/t-dependence with the step at A/t = 2to- 
Next, we consider the case of a highly asymmetric cou- 
pling to the leads, (77 + 1/77) ^> 1. For Ap > 2t , there is 
a competition between two types of processes of inelastic 
cotunneling. One is the thermal equilibration of the DD, 
due to inelastic cotunneling into the same lead; and the 
second one is the heating effect of the DD, due to inelastic 
cotunneling from the left lead to the right lead, provided 
PL > PR- The strength of the former effect is propor- 
tional to exp(-2£ /T), for T < 2t . The latter effect is 
proportional to (T/2t )rj/(l + rj 2 ), for Ap = 2t . As a 
function of T, the crossover occurs at the energy scale 
Xh = 2to/w(r] + I/77), where w(x) — \m{x ln(x... ln(x))) 
for x > e. At T 3> Th, the DD is in thermal equilibrium 
with 7 = p-/p+ = exp(-2t /T). At T ~ T h , the heating 
in the DD is governed by both types of processes, and the 
ratio 7 depends on both T and A/i, as given by Eq. 1)620. 
At T <C Th, we are in the strong heating regime, domi- 
nated by processes of inelastic cotunneling from one lead 
to the other. Here, we have 



( v T 

_ J 2t (l +I ?) 2 ' 

T — S A A1 -2t n 

^ A/i+2t (l+ t; +l/77) ' 



|A M -2< | «T, 
A/x - 2i > T. 



(66) 



In this regime, we can extract the asymmetry parameter 
77, for A/i — 2to 3> T, in the following way 



4?/ 



(l + r,)2 



1- A 



A/j-2t 
4tn 



.4 



2i 

AGdAp'' 



(67) 
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FIG. 9: (a) Differential conductance G vs bias in the cotun- 
neling regime for the N — 2 CB valley. The dashed curve 
corresponds to the asymptotic value of G at low tempera- 
tures. The solid curves are for a finite T = (0.2, l,2)Th, 
where Th is the characteristic temperature of the strong heat- 
ing regime, see text. For the calculation we used: <S = 0.6, 
cj> = 0.316 (T h ~ 0.23J), 77 = 1, E-{2) = E+{2) = E c , and 
Tl^/^c = 0.1. (6) The quantity A/[l ~ A(Afi ~ J)/2J], 
with A = - (J/ AG) (dG/dAu) , where AG = G(Afx) - G(oo), 
is plotted for the three solid curves of G(A/i) shown in Fig.|5Jx. 
The solid curve, corresponding to T = 0.2Th, saturates with 
good precision at the value of 8/(2 + ft) = 2(1 — cj>) 2 /(l + (j) 2 ), 
whereas the other solid curves (not in the strong heating 
regime) have a smaller saturation value. The dotted curves A, 
B, and C, are plotted for the T = 0.2Th curve and correspond 
to choosing the value of G(oo) as shown in Fig.|5Jj by dotted 
horizontal segments. Only for a correct choice of G(oo) the 
curve has a plateau and saturates at a finite value. 



where AG = G(A/i) - G(oo). The value of G(oo) is the 
value of G at A/i - 2t > 2£ (1 + ■q) 2 jr). We note that 
Eq. (|67(l holds also for 77 ~ 1, then the energy scale Th 
coincides with 2io- 



V. COTUNNELING IN THE N = 2 COULOMB 
BLOCKADE VALLEY 



tracing out the spin degeneracy of the triplet the co- 
tunneling rates are given by Eq. (JHTJl with M l n l m given 111 
Appendix [E] 

The following discussion is similar to the one in Sec- 
tion llVI The heating effect of the DD is described by the 
ratio 



The width of the N = 2 CB valley is of the order of U H - 
The energy scale of interest here is the exchange J <C 
Uh- Similarly to Section lTVI we calculate the cotunneling 
rates with N — 2 electrons in the DD, using the formula 
with Tin given by 



-in 



E-{2) 



' E + (2) 



, Wal Cl h a -(68) 



Here, we assumed 2t$ -C E_(2) and the energy split- 
ting in the N = 3 sector to be much smaller than E+(2) 
as well as A/j,T < E±{2). The DD states \m), \n) in 
(|59|l are now the singlet-triplet states @ and the cor- 
responding energies are taken from JSJ and 10. After 



WTS _ 
WST 

9(- J + Am) + 8(— J - A/j) + k0(- J) 

e(j + A/j) + e(j- A/j) + K e(j) ' 1 ' 



where k = (rj + l/ij)(l + cf>) 2 / (1- <j)) 2 . The DD population 
probabilities are given by ps = 1/(1-1-/3) and /?t = /?/(! + 
f3). The formulas l|63|l . (|64|l and (|65|l apply to this case 
also, provided we substitute the indices + by S and — by 
T. The differential conductance G, at T <C J, has a step- 
like A/i-dependence, with the step occurring at A/i = J. 
We plot G vs Am for different temperatures in Fig.|5Ji. 

The strong heating regime was considered previously^!" 
for the case 77 = 1. We present the results for arbitrary 
77 here. The energy scale of the strong heating regime for 
k 3> 1 is given by Th = J/w(k), where the function w(x) 
was introduced in Section IIVI The asymptotes of (3 in 



11 



the strong heating regime, T -C Th, are as follows 



3T 



= 



J(2+k) ' 

3(A^-J) 
A^i+,7(1+k) 



|A/i - J| « T, 
A/i - J > T. 



(70) 



For A/i - J » T, the following equality holds in the 
strong heating regime, 



.4 



2 + k l-A^f^' 



.4 



J rfG 
AGdA/? 



(71) 



where AG = G(A/i) — G(oo), with G(oo) given by G at 
A/i — J ^> (2 + k) J. Eq. (f7Tj) allows us to extract k, and 
hence, to extract the interaction parameter (f>, provided 
the asymmetry parameter n is known (e.g. from a similar 
procedure as explained in Section llVjl . We illustrate how 
this procedure works in Fig. The r.h.s. of Eq. I|71|l 
is plotted for the curves in Fig. [5Ji. With lowering the 
temperature, the plateau value saturates at 8/(2 + k), as 
one enters the strong heating regime; compare the three 
solid curves in Fig.[2>. If the value of G(oo) is not known, 
then one can proceed in the following way. Starting with 
a lower bound of G(oo), which can be e.g. G(oo) = 0, 
a set of curves is plotted for values of G(oo) increased 
each time by an offset value. The curves can be divided 
into two classes: (i) curves which have no divergence for 
A/i — J > T, have a maximum, and saturate at zero; and 
(ii) curves which monotonically increase, or even diverge, 
for biases A/i — J > T in the available (measured) range 
of A/i. The separatrix of these two classes has a mono- 
tonic dependence with a plateau at a non-zero value, and 
it corresponds to the asymptotic value of the cotunnel- 
ing conductance G(oo). The dotted curves in Fig. |5f), 
denoted as A, B, and C, are plotted for values which 
might mistakenly be assigned to G(oo), and are related 
to the curve at T = 0.2Th. The curves A and B belong 
to the class (i), and the curve C to the class (ii); the 
separatrix is the solid curve at T = 0.2Th. The values 
of G(oo) taken for the curves A, B, and C, are shown in 
Fig. Eh by dotted lines, whereas the true value of G(oo) 
by a solid line. Next, assuming that the measurement of 
G(A/i) has an error bar, we note that the value at the 
maximum in curves of class (i) provides a lower bound 
for 8/(2 + k). 

Finally, we note that the same physics holds true for 
the case J < (triplet ground state), which can be re- 
alized by applying a magnetic field perpendicular to the 
2DEG plane. Here, one should replace J — ► |J| and 
-> 9//3 in (25, and 8/(2 + k) -> 8/(6 + 3k) in the l.h.s. 
of i|71|) . Eq. I|69(l remains valid for this case. 



VI. COTUNNELING-ASSISTED SEQUENTIAL 
TUNNELING 

In this section we consider the interplay of sequential 
tunneling and cotunncling close to the sequential tun- 
neling peaks at finite bias. We find that at low tem- 
peratures, T <C T , the heating effect due to cotun- 
neling provides population to the excited states, from 



which a subsequent sequential tunneling can occur. Such 
a cotunneling- assisted sequential tunneling produces new 
features (peaks/dips) in G versus A/i and SE. The en- 
ergy scale To is given by To = 2to/ ln(2£o/r.R) for the 
N = 1 CB side, and T = J/ln(J/r L ) for the N = 2 
CB side; here, we assumed symmetric biasing of the DD 
with A/i£ = A/i/2 > 0. For G versus A/i we find a peak 
at A/i/2 = \5E\ - 2t > 2t on the N = 1 CB valley side, 
as well as a peak/dip at A/i/2 = \SE\ — 2t + J. On the 
N = 2 CB valley side, we find a peak in G versus A/i at 
A/i/2 = \SE\ - J > J. 

We proceed with considering the N = 1 CB valley side 
and for the sake of simplicity of the following expres- 
sions we assume J -C 2i <C U12, which corresponds to 
weakly coupled dots. We consider a position in the CB 
valley close to the N = 1,2 sequential tunneling peak, 
T_(l) ^> E+(l), but still far enough to be able to ap- 
ply a bias A/i > 2to and to have £+(1) — A/i^ > 2to. 
Since A/i^ is comparable to — SE = E + (l), the cotun- 
neling rates obtained in Section IIVI are not valid here 
and should be modified as to account for the energy de- 
pendence of the tunneling density of states in the bias 
window. We replace Eq. (|6Ufl by 



Tin = - 



E 



C l'k'a' C lk<Ji 



(72) 



where U{ = E{2) - E(l) - m, and A was defined below 
equation (JBDJ. Formula J72J) is valid for T <C U l + . Next, 
we calculate the cotunneling rates using the golden rule 
expression (|59|l and trace out the spin degeneracies ac- 
cording to Eq. 122fl . We obtain the cotunneling rates 
close to the sequential tunneling peak, 



1 



w nm (l ,1) = —v [ W - Ul+E 



M 



(73) 



for E mn + A/i/;/ > 0, and w nm (l',l) suppressed by 
ex.p[(E mn + Afiw)/T], for E mn + A/i«/ < 0. Here, we 
have assumed T <§; 2t a and \E mn + A/i«/| >• T. The 
quantities M l n l m are obtained from M l n l m in Appendix iDl 
by setting U_ — > 00 and omitting the l/iT 2 . denomi- 
nators. We note that the diagonal rates w nn (l,l) do 
not enter our further calculation, and hence the case 
E mn + Afiw = refers only to the vicinity of A/i = 2io, 
where the value of the rate is proportional to T. For this 
case (\E mn + Afi w \ <C U l + ), one can use the expressions 
derived in Section El setting £L(1) > E + (l). 

Next we solve the master equation, similar to Sec- 
tion ^^^J including also the cotunneling rates. We find 
that, on the N — 1 CB valley side, it suffices to take into 
account only the rates 1(751 . i.e. between |+) and |— ), and 
neglect the cotunneling between \S) and |T) (they give 
higher order correction). From the competition of the 

cotunneling rate w 1- and the sequential tunneling rate 

Ws,+ we deduce the energy scale T = 2i / ln(2io/fifj), 
valid for A/i > 02£. For temperatures T <C Tq, the ratio 
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FIG. 10: Differential conductance G versus bias for: (a) SE = 
— 1 and (b) SE — 1. Here, we use the following parameters: 
t = 0.1, J = 0.025, S = 0.2, 4> = 0.6, = Tr — 0.001, 
T — 0.001, A/j,l,r = Ap/2. The conductance is scaled by 
Gr = e 2 r R /hT on both plots. 



7 = P-/P+ is given by 

W+, S W S , 



7 = W-+ u> + - + 



+ 



W+.s + W-,s W+, T + W-, T 



(74) 

This straightforwardly yields a non- vanishing population 
in the N = 2 sector, 



: 7 



W s ,- 



+ 



It 



r,- 



v W+,s + VF- S ' W +I T + W-,T , 

with the balance between |T) and jS 1 ) given by 



W s ,- (W+, T + W-, T ) 



(75) 



(76) 



The physical interpretation of equation (|74|l is as follows. 
In the bias range given by Ui > 2to, the occupation 
of the DD is determined by the cotunneling processes, 

with 7 = w |_/k; + _, see Section llVl For biases A/t > 

2to, the population p_ increases with increasing the bias, 
see l|66f) . and can reach a value comparable to p + . A 
further increase of the bias brings us to the vicinity of 
the point ?7f = 2io- Here, a new channel of relaxation 
opens, namely |— ) — > \S) — ► |+) (and also |— ) — > |T) — > 
|+) at = 2to — </)• Since the reverse sequence is 
forbidden due to energy conservation, the level |— ) is 
emptied efficiently and one has w 1. A small non- 
equilibrium population persists in the N — 2 sector and 
in the level |— ). This allows for a sequential tunneling 
current from the excited DD states. 

The current through the DD consists of two parts, 



i = / scq + r 



(77) 



where / scq is the sequential-tunneling part given by 
Eq. H48fl with ps,T calculated above. The cotunneling 
part 7 cot can be calculated with the Eq. I|b3[l and the 
rates f75|). Both terms, 7 scq and 7 cot in l(77|l. contribute 
with the same order of magnitude to the current and 



differential conductance. In Fig. 1 10b . we plot the bias 
dependence of the differential conductance G on the left- 
hand-side of the sequential-tunneling peak. The effect 
discussed above results in two features, namely (i) a peak 
at Ap/2 = \SE\ — 2t and (ii) a peak/dip at A/i/2 = 
\SE\ — 2t + J (both positions are given for T — 0). The 
latter peak goes into a dip with increasing the inter-dot 
coupling 30 . The positions of the peaks are temperature 
dependent. With increasing T both peaks shift to the 
left; the peak (i) shifts by 5fi cx Tln[2t /ut 2 R (l + l/2r?)] 
and it has a width of the same order of magnitude, the 
peak (ii) shifts by <5/i ~ T and has width ~ T. This scal- 
ing behavior with varying = 7rt4|j comes about from 

the competition between the cotunneling rate and 

sequential tunneling rate Ws,- in equation l|74|) : for the 
peak (ii) it is not present. We note that the peak (ii) is 
solely due to dI seq /dAfi, whereas the peak (i) is present 
in both dI scq /dAii and d/ cot /rfA/i. 

Next, we proceed with considering the N = 2 CB valley 
side close to the N = 1,2 peak (£L(2) < £+(2)). We 
assume J <C 2i <C Uh, which is usually the case for 
double dots, and consider a position in the valley such 
that one can satisfy A/i > J and E_(2) — A/i^ > J. 
We follow a derivation close to the case discussed above. 
Instead of Eq. H72|) we use 



T 



VI 



- - ^2 **n*i'«' 



(78) 



where U l _ = E(l) ~ E{2) + m and A = J2 a l+> CT > (+> 
To calculate the golden rule rates we use Eq. I|73|l with 
U_ — * U l + and Ai l n l m given in Appendix El From 
solving the master equation we find results similar to 
Eqs. (|74|1 — (|76|l . Namely, we find that in the tempera- 
ture regime T < T , with T = J/ ln( J/vt\) for Afi > 
and To = J / h\(J / vt 2 R ) for A/i < 0, the cotunneling rate 
wts provides non-equilibrium population to the state |+) 
for U* < J, where U R = E_(2) - Ap R . The ratio 
[3 = pt/ps in this regime is as follows 



[3 = W T S W S T 



Ws, + W +i t 



The population in the N 



1 



W s ,+ + W T ,+ 
• 1 sector is determined by 

W +<T 



(79) 



W s , + + W t ,a 



(80) 



and belongs to the |+) state, i.e. 7 = 0. The cur- 
rent though the DD is given by Eq. I|77|) . The sequen- 
tial tunneling part J seq is given by Eq. (|49|l . where one 
should assume p_ = 0, and />+ = 1/t <^ 1 given 
by Eq. (|80|l . The cotunneling part I cot is given by 
Eq. (|rJ5|> . where one should use occupation probabilities 
determined by Eqs. I|79l) and l|80() . and cotunneling rates 
as discussed after Eq. l(7S|) . In Fig. 110b. we plot G vs 
A/x/2 = A/.ir > on the left-hand-side of the sequential- 
tunneling peak. The new feature is a peak, which occurs 
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at A/z/2 = SE — J for T = 0. With increasing T, it moves 
to the left by an amount <5/x oc T\v\J jvt\{\ + rj/2)], and 
has a width of the same order of magnitude. We note 
that the behavior of this peak is analogous to that of the 
peak (i) in the previously discussed case (see above). 

To summarize, we have analyzed an effect originating 
from an interplay of sequential tunneling and (inelastic) 
cotunneling, and have specified the regimes when addi- 
tional features in transport through a DD can be ob- 
served. 



VII. CONCLUSIONS 

We have analyzed theoretically the transport spec- 
troscopy of a symmetric DD attached to leads in series. 
Our motivation was to find ways of characterizing the 
DD in dc transport measurements. For this, we have 
recast the main results of the Hund-Mulliken method 19 , 
introducing a description of the DD in terms of a set 
of parameters, {io, J, 4>, S, U12, Uh, So>o}i which can be 
referred to in experiments as phenomenological param- 
eters, see Sec. |H] Direct access to these parameters is 
desirable for building spin-based qubits using quantum 
dots, see Ref. |3ll Using a master equation approach, 
we have described the transport and the non-equilibrium 
probability distribution in the DD both in the sequential 
tunneling and cotunneling regimes. We have specified a 
number of "non-universal" regimes, which reveal infor- 
mation about the DD parameters. We summarize our 
main results below. 

In the sequential tunneling regime, see Sec. lIIII the dif- 
ferential conductance at a finite bias G(A/j,) has satellite 
peaks of sequential tunneling with respect to the main 
peaks of the CB diamond, see Figs.[7||Hl and[S] The ex- 
change coupling constant J appears as a peak separation 
in G(Afx) both on the N = 1 and N = 2 CB sides. The 
tunnel splitting 2tg can be extracted from G(Afi) on the 
N = 2 CB side, see Fig. EJ We note that if the two dots 
are detuned by some energy AE < 2to, then one can re- 
place 2to — > \J 4£q + AE 2 for the transport spectroscopy 
oftheDD. 

In the cotunneling regime, see Sees. M\ and E\ the 
exchange coupling constant J can be observed on the 
N = 2 CB side as the bias value at which the inelastic 
cotunneling turns on. A step in G(A/x) occurs at \Afj,\ = 
\J\. Similarly the tunnel splitting 2io can be extracted 
from the N = 1 CB side. In Sec. II VI we have introduced 
a strong heating regime where the bias dependence of 
G(A/i) allows one to extract the asymmetry parameter 
V = |*fl| 2 /l*i| 2 on the N = 1 CB side. An additional 
relation involving r/ and 4> can be obtained on the N = 2 
CB side, see Sec. H3 

As an alternative to transport measurements we have 
considered the use of a charge detector (QPC) close to 
one of the dots, see Sec. IIII CI We found that measur- 
ing the average charge on the DD allows one to extract 
J and 2£q, similarly to the conductance measurements 



in the sequential tunneling regime. Moreover, additional 
relations between the DD parameters can be obtained 
with sensitive QPC measurements. We have also consid- 
ered the back action of the QPC onto the DD and found 
that the QPC induces relaxation of the DD states with 1 
electron. We accounted for this relaxation in the master 
equations in Sec. IIII CI 

Finally, we have analyzed a combined mechanism 
of sequential tunneling and cotunneling, see Sec. IVII 
We found that inelastic cotunneling can provide non- 
equilibrium population probability to the excited DD 
states, which can then allow for sequential tunneling via 
an excited DD state. Accounting for this effect results in 
additional satellite peaks in G(Afi) at low temperatures. 



APPENDIX A 



Using the definition 



V2 



nd 



(Al) 



we rewrite (0 in terms of the new operators di a as follows 



(i-^)(4 T 4i+4r4i)}i ) . 

\T+) = 4 T 4 t |0) , |T_> = d[/ Rl |0) , (A2) 



1 

3 V 



Similarly, for (|10fl and l|ll|) we have 



(A3) 



= ■ " 1)( ^t4i - 4A)+ 



2v^T^ 

)} |0> • (A4) 
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APPENDIX B 



The sequential tunneling rates W l Mm calculated ac- 
cording to the formula l|19(l are: 



w U + ,d 



W \s\-A) 



f( E \s) - E \+,T) ~ Mi) 



27T 4> 2 \tL~\ 2 



w (s\-,d = t u i 



2 m^-E^-u), 

27T <^_| 2 



W, 



W (T + \ + ,l) 



2n 

Ti 
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^-v\t L _\ 2 f(E\ T+) -E M) 



W\ T+ 


1- T) 
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| 2 /(^|T+> 


~ E \--1) 


- w) . 


W (T + 
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= o, 








w[ T _ 


l+.T) 


= o, 








W\ T _ 
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2tt . 

= ¥^- 


\ 2 f{E\T_) 


~ E \ + ,i) 


- Mi) > 


W\ T _ 


l-.T) 


-o, 








w\ T _ 




2tt . 
= T ^ + 


| 2 /(^|T_) 


-E\- tl ) 


- Mi) , 




+,T> 


2tt |t,,_ 
ft 17 2 


I 2 


~ E \+A) 


- Mi) , 


W\ Ta \ 


+,1) 


2tt |tj,_ 
ft ^ 2 


I 2 


-%,!> 


- Mi) j 


W U\ 


- T> 


_ 2tt |t i+ 
ft ^ 2 


| 2 

f( E \T ) 


~ ^1- T) ' 


- Mi) > 


W\t \ 


- i> 


27T jt £+ | 2 

" ft " 2 /(i? l To > 


-^1-4) 


- Mi) > 



where f(E) = 1/ [1 + exp(E/k B T)]. The rates for the 
reverse transitions can be obtained from the above for- 
mulas by replacing f(E) — ► 1 — f(E), satisfying (|2*T|l . 



APPENDIX C 



Solving (d) and 123, we find for {25, {23, and QEJ 

the following expressions: 



(Wt,-W-. s + W-^rgV) (Wg,+ + W T ..+) + (V^g, + W+,T + Wr,+W+,s) (W T - + Wg ,-) 
(W +>S W S _ + gWg,+) (W+, T + W_, T ) + (W+,tWt,- + W_, t Wt,+) (W+,g + W_ <?) 
W T , + V^ + ,g (Wg - + W T ,-) + W T -W-, S {W s ,+ + W T ,+ ) 
W+,tW s ,+ {Ws- + W T -) + W-, T W S - (W s ,+ + W T ,+) ' 
W_ gWg,+ (W+,t + W-, T ) + W-, T W T ,+ (W+,g + VK_,g) 
W+, S W S ,_ (W+, T + VF_,t) + W + , T W T ,- (W+,s + W_ <?) ' 



(CI) 
(C2) 
(C3) 



(Wg+ - W S -)(W+t + W-r) + (Wr+ - W T -){W +S + W_g) 
W+gWg_(VK fT + W_ T ) + W +t Wt-{W +s + W_ s ) 



(C4) 



15 



APPENDIX D 
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3 ^ 










-U+ 


+ 2Ul; 



2 + 3|tL 2 ^r |2 - j^jfi {t* R+ t R -t L+ t* L _) , (di) 



^ = ( TSr + T^TrV + 7TT~l2V2 772" ) 2 l^*«-| 2 + - ^« fet^O , (D2) 

where f/_ = £L(1) and U+ = E + (l), and 3? is real part. For change U+ «=* in Eq. |jD3j); for M^m> set 

L — > i?. 
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where U- = £L(2) and [7+ = £+(2). 
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